##Opis danych i cel projektu Dane przedstawiają inforamcje o pacjentach z niewydolnością serca.Pacjenici zostali podzieleni na 4 kategorie: 1)Palacz, niewydolność śmiertelna 2)Palacz, niewydolność nie śmiertelna 3)Nie palacz, niewydolność śmiertelna 4)Nie palacz, niewydolność nie śmiertelna

Zmienne v1,..,v7 określają: 1)wiek 2)poziom kinazy kreatynowej we krwi (mcg/l) 3)procentowy poziom krwi opuszczającej serce przy każdym skurczu 4)poziom płytek krwi 5)pozim kreatyniny w surowicy (mg/dl) 6)poziom sodu w surowicy (mEq/l) 7)okres kontrolny (w dniach)

Źródło:https://doi.org/10.1186/s12911-020-1023-5 https://archive.ics.uci.edu/dataset/519/heart+failure+clinical+records

##Analiza danych jednowymiarowych

head(dane)
##   Kategoria V1   V2 V3     V4 V5  V6 V7
## 1         3 75  582 20 265000  2 130  4
## 2         3 55 7861 38 263358  1 136  6
## 3         1 65  146 20 162000  1 129  7
## 4         3 50  111 20 210000  2 137  7
## 5         3 65  160 20 327000  3 116  8
## 6         1 90   47 40 204000  2 132  8
any(is.na(dane))
## [1] FALSE

Na początku zaprezentuję hiostogramy oraz podstawowe statystyki opisowe dla wszytkich zmiennych

"podsumowanie"
## [1] "podsumowanie"
apply(dane[,-1],2,summary)
##               V1        V2       V3     V4      V5       V6       V7
## Min.    40.00000   23.0000 14.00000  25100 1.00000 113.0000   4.0000
## 1st Qu. 51.00000  116.5000 30.00000 212500 1.00000 134.0000  73.0000
## Median  60.00000  250.0000 38.00000 262000 1.00000 137.0000 115.0000
## Mean    60.83612  581.8395 38.08361 263358 1.41806 136.6254 130.2609
## 3rd Qu. 70.00000  582.0000 45.00000 303500 1.00000 140.0000 203.0000
## Max.    95.00000 7861.0000 80.00000 850000 9.00000 148.0000 285.0000
"odchylenie standardowe"
## [1] "odchylenie standardowe"
sort(apply(dane[,-1],2,sd), decreasing = TRUE)
##           V4           V2           V7           V1           V3           V6 
## 97804.236869   970.287881    77.614208    11.894809    11.834841     4.412477 
##           V5 
##     1.043906
"wariancja"
## [1] "wariancja"
sort(apply(dane[,-1],2,var), decreasing = TRUE)
##           V4           V2           V7           V1           V3           V6 
## 9.565669e+09 9.414586e+05 6.023965e+03 1.414865e+02 1.400635e+02 1.946996e+01 
##           V5 
## 1.089740e+00
"skośność"
## [1] "skośność"
sort(apply(dane[,-1],2,skewness), decreasing = TRUE)
##         V2         V5         V4         V3         V1         V7         V6 
##  4.4406886  4.2072326  1.4549746  0.5525927  0.4203739  0.1271606 -1.0428705

Największe wariancje mają zmienne v4 oraz v2, zmienne te mają też największe średnie oraz odcylenie standardowe. Wszystkie zmienne poza zmienną v6 mają dodatnią skośność.

for(i in 1:7) hist(dane[,i+1],main=colnames(dane)[i+1])

##Rozkład wielomodalny i obserwacje odstające

Rozkład normalny przypominają zmienne V6 oraz V4. Natomiast jeżeli chodzi o wielomodalnośc to przetestuje zmienne V3,V7,V2 oraz V5.

shapiro.test(dane$V6)
## 
##  Shapiro-Wilk normality test
## 
## data:  dane$V6
## W = 0.93903, p-value = 9.215e-10
shapiro.test(dane$V4)
## 
##  Shapiro-Wilk normality test
## 
## data:  dane$V4
## W = 0.91151, p-value = 2.883e-12

Żadna zmienna nie ma rokładu normalnego

silverman.test(dane$V3,k=1)
## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.04004004
silverman.test(dane$V7,k=1)
## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.001001001
silverman.test(dane$V2,k=1)
## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.07007007
silverman.test(dane$V5,k=1)
## Silvermantest: Testing the null hypothesis that the number of modes is <=  1 
## The resulting p-value is  0.04204204

p-value poniżej 5% mają zmienne V3,V7,V5. Oznacza to, że mają roakłady wielomodalne. Za pomocą silvermanplot sprawdzę ilo modalne są te zmienne

silverman.plot(dane$V3)

silverman.plot(dane$V7)

silverman.plot(dane$V5)

Zwykresu zmienne V3 oraz V5 są trójmodalne natomiast zmienna V7 jest cztero lub pięcio modalna.

Teraz korzystając z testu Grubbsa sprawdzimy obserwacje odstające.

apply(dane[,-1],2,function(x) grubbs.test(x))
## $V1
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 2.87217, U = 0.97222, p-value = 0.576
## alternative hypothesis: highest value 95 is an outlier
## 
## 
## $V2
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 7.5021, U = 0.8105, p-value = 4.647e-13
## alternative hypothesis: highest value 7861 is an outlier
## 
## 
## $V3
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 3.54178, U = 0.95776, p-value = 0.05195
## alternative hypothesis: highest value 80 is an outlier
## 
## 
## $V4
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 5.99812, U = 0.87887, p-value = 9.134e-08
## alternative hypothesis: highest value 850000 is an outlier
## 
## 
## $V5
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 7.26305, U = 0.82239, p-value = 3.95e-12
## alternative hypothesis: highest value 9 is an outlier
## 
## 
## $V6
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 5.35423, U = 0.90348, p-value = 6.145e-06
## alternative hypothesis: lowest value 113 is an outlier
## 
## 
## $V7
## 
##  Grubbs test for one outlier
## 
## data:  x
## G = 1.99370, U = 0.98662, p-value = 1
## alternative hypothesis: highest value 285 is an outlier

Zmienne V1,V3,V7 nie mają obserwacji odstających

##Macierz korelacj, redukcja wymiaru

Macierz Pearsona

heatmaply_cor(cor(dane[,-1]))

Ogólnie dane są ze sobą słabo skorelowane. W delikatnym stopniu można zobaczyć korelację danych V1 i V5, oraz V6iV3.

heatmaply_cor(cor(dane[,-1],method = 'kendall'))

Macierz korelacji Kendalla jest bardzo podobna do macierzy Pearsona. Można na niej zobaczyć jednak odrobinę słabsze skorelowanie danych poza danymiV1 i v5 które są ze sobą skorelowane minimalnie bardziej.

prcomp(dane[,-1])->pca.dane
summary(pca.dane)
## Importance of components:
##                              PC1     PC2   PC3  PC4   PC5   PC6   PC7
## Standard deviation     9.780e+04 9.7e+02 77.66 12.1 11.23 4.305 1.009
## Proportion of Variance 9.999e-01 1.0e-04  0.00  0.0  0.00 0.000 0.000
## Cumulative Proportion  9.999e-01 1.0e+00  1.00  1.0  1.00 1.000 1.000
df.pca<-data.frame(pc1=pca.dane$x[,1],pc2=pca.dane$x[,2],pc3=pca.dane$x[,3],kl=as.factor(dane$Kategoria))
plot(x=df.pca[,1],y=df.pca[,2], col=as.factor(dane[,1]))

plot_ly(df.pca,x=~pc1,y=~pc2,z=~pc3,color=~kl,type='scatter3d')
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
sort(pca.dane$rotation[,1])
##            V1            V5            V6            V7            V3 
## -6.359893e-06 -4.606192e-07  2.802784e-06  8.343479e-06  8.733849e-06 
##            V2            V4 
##  2.427180e-04  1.000000e+00
sort(pca.dane$rotation[,2])
##            V2            V6            V5            V4            V3 
## -9.999990e-01 -2.640333e-04  2.199171e-05  2.427134e-04  5.596420e-04 
##            V7            V1 
##  7.733495e-04  9.861072e-04
sort(pca.dane$rotation[,3])
##            V7            V3            V6            V2            V4 
## -9.993492e-01 -6.271944e-03 -5.006069e-03 -7.403701e-04  8.810713e-06 
##            V5            V1 
##  1.866727e-03  3.510970e-02

Pierwsza składniowa zawiera prawie 100% wariancji.PC2 jest zdominowane przez zmienną V2, PC3 przez zmienną V7.

##T-sne

dane.tsne<-Rtsne(dane[-1],theta=0,dims=3,perplexity=5)
dane.tsne.2<-Rtsne(dane[,-1],theta=0.5,dims=3,perplexity=20)
summary(dane.tsne$Y)
##        V1                V2                V3         
##  Min.   :-58.489   Min.   :-48.003   Min.   :-48.013  
##  1st Qu.: -8.766   1st Qu.:-22.906   1st Qu.:-17.219  
##  Median : -1.094   Median :  3.105   Median : -1.518  
##  Mean   :  0.000   Mean   :  0.000   Mean   :  0.000  
##  3rd Qu.: 15.767   3rd Qu.: 12.119   3rd Qu.: 17.640  
##  Max.   : 64.226   Max.   : 53.072   Max.   : 54.333
summary(dane.tsne.2$Y)
##        V1                 V2                 V3         
##  Min.   :-11.9476   Min.   :-18.7120   Min.   :-13.299  
##  1st Qu.: -6.8698   1st Qu.:-11.6139   1st Qu.: -9.430  
##  Median :  0.1425   Median : -0.5326   Median : -2.878  
##  Mean   :  0.0000   Mean   :  0.0000   Mean   :  0.000  
##  3rd Qu.:  6.9093   3rd Qu.: 10.4147   3rd Qu.:  9.972  
##  Max.   : 10.4579   Max.   : 21.4740   Max.   : 19.852
dane.tsne.df<-as.data.frame(dane.tsne$Y)
dane.tsne.2.df<-as.data.frame(dane.tsne.2$Y)
plot_ly(dane.tsne.df, x = ~V1, y = ~V2, z = ~V3,color = dane$Kategoria, type = 'scatter3d')
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
plot_ly(dane.tsne.2.df, x = ~V1, y = ~V2, z = ~V3,color = dane$Kategoria, type = 'scatter3d')
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Dane są ułożone chatycznie. Ze wzrostem parametru perplexity zbijają się w grupy, a poszególne grupy oddalają się od siebie.

##Analiza skupień Zacznę od ilości klasrów równej ilości kategorii

kmeans(dane[,-1],centers=4)->km.dane.2
df.pca<-data.frame(pc1=pca.dane$x[,1],pc2=pca.dane$x[,2],pc3=pca.dane$x[,3],kl=as.factor(km.dane.2$cluster))
plot_ly(df.pca,x=~pc1,y=~pc2,z=~pc3,color=~kl,type='scatter3d')
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Wygląda na to, że k-średnich ładnie podzieliło wykres, ale upewnimy się jeszcze

table(dane$Kategoria, km.dane.2$cluster)
##    
##      1  2  3  4
##   1  5  6  3 16
##   2 17 16  3 30
##   3 19 17  2 28
##   4 20 31  6 80

mogło by być lepiej

wss <- rep(NA,9)
for(i in 1:9) wss[i] <- kmeans(dane[,-1], centers=i+1)$tot.withinss
plot(1:9,wss)

Z wykresu wynika, że 3 to optymalna liczba klastrów. Przeprowadzę więc klastrowanie dla tej liczby

kmeans(dane[,-1],centers=3)->km.dane.3
df.pca.3<-data.frame(pc1=pca.dane$x[,1],pc2=pca.dane$x[,2],pc3=pca.dane$x[,3],kl=as.factor(km.dane.3$cluster))
plot_ly(df.pca.3,x=~pc1,y=~pc2,z=~pc3,color=~kl,type='scatter3d',mode='markers')

Znowu wygląda ładnie, ale sprawdzimy to

table(dane$Kategoria, km.dane.3$cluster)
##    
##      1  2  3
##   1 10 16  4
##   2 27 32  7
##   3 30 29  7
##   4 44 80 13

Znowu nie wygląda to dobrze, najwięcej wyników wrzuciło do klastra 3

Przeprowadze teraz klastrowanie aglomeracyjne oraz podziałowe.

plot(agnes(dist(dane[,-1], method='minkowski',p=1),diss=T))

plot(diana(dane[,-1]))

Wyniki są zbliżone ##Klasyfikacja pod nadzorem- lasy losowe

rfcv(trainx = dane[,-1],trainy = as.factor(dane$Kategoria))->rf.dane
rf.dane
## $n.var
## [1] 7 4 1
## 
## $error.cv
##         7         4         1 
## 0.4916388 0.5050167 0.5451505 
## 
## $predicted
## $predicted$`7`
##   [1] 3 3 3 3 3 3 3 3 3 1 3 3 3 3 4 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 3 3 3 3
##  [38] 3 3 3 3 3 3 3 1 3 1 3 3 4 1 1 3 4 3 3 3 3 3 1 3 3 3 3 3 3 1 3 3 3 3 2 2 2
##  [75] 1 3 4 4 4 2 4 4 3 3 4 4 4 4 4 2 4 4 4 3 4 4 4 2 4 4 4 4 4 4 4 4 4 2 4 4 4
## [112] 4 3 4 4 4 4 3 4 3 3 4 4 4 3 4 3 4 2 4 4 3 4 4 4 4 4 3 4 2 4 2 4 4 4 4 4 4
## [149] 3 4 4 4 4 4 4 3 2 4 3 4 4 4 4 2 4 4 4 3 4 4 4 2 2 2 4 4 4 4 2 4 4 1 4 3 4
## [186] 3 4 1 4 4 3 4 4 4 3 4 4 2 4 4 4 4 4 3 4 4 4 3 4 4 4 4 4 4 4 4 4 2 4 2 4 4
## [223] 4 4 4 4 4 2 3 4 4 4 4 4 4 2 4 2 4 4 4 4 4 4 2 4 4 4 4 4 3 4 4 4 2 4 4 4 4
## [260] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 3 3 4 4 4 4 4 4 4 4 4 4 4 4 2
## [297] 2 4 4
## Levels: 1 2 3 4
## 
## $predicted$`4`
##   [1] 1 1 3 3 3 3 1 3 1 1 3 3 4 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 3 3 3 3 3
##  [38] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 3 3 2 3 3 1 3 3 3 2 2 3 1 3 1 3 2 2 2 2
##  [75] 1 3 4 4 2 2 2 4 3 3 4 4 4 4 4 3 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 3 4 2 4 4 4
## [112] 4 3 4 4 4 4 3 4 4 4 4 4 2 4 4 3 4 4 2 4 3 4 4 3 4 4 3 4 2 3 2 4 4 3 4 4 4
## [149] 3 4 4 2 4 4 4 1 4 1 4 4 4 4 4 2 2 3 4 3 2 4 4 2 4 3 4 4 4 3 2 4 2 1 4 3 1
## [186] 3 4 1 4 4 3 2 4 4 3 4 4 2 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 3 4 4 4 2 4 4 3 4
## [223] 4 4 3 4 4 4 3 2 4 4 2 4 4 2 2 2 4 4 2 2 2 4 4 4 4 4 4 4 2 2 4 4 2 4 4 4 4
## [260] 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 2 4 4 4 2 4 4 4 4 2 4 4 2 2
## [297] 2 2 2
## Levels: 1 2 3 4
## 
## $predicted$`1`
##   [1] 1 1 3 1 1 3 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 1 3 1 4 4 1 3 1 3 1 3 3 3 4 4 4
##  [38] 4 3 3 3 1 3 3 1 1 1 3 1 3 3 3 3 3 3 2 2 2 2 3 3 2 3 3 1 4 1 3 3 3 3 1 2 1
##  [75] 1 1 4 2 2 2 2 4 2 3 4 4 4 4 4 4 4 2 3 4 4 4 4 4 4 4 4 2 4 4 2 2 1 2 1 4 2
## [112] 3 1 2 2 2 4 4 4 4 4 4 4 4 4 3 4 2 4 2 2 2 4 4 4 2 4 4 4 4 4 4 4 4 4 2 2 3
## [149] 4 3 4 2 4 2 4 4 4 4 4 4 2 4 2 2 3 3 3 2 3 2 4 2 2 2 2 2 4 2 2 4 4 4 1 1 3
## [186] 3 3 3 3 4 4 4 4 3 2 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 3 4 3 4 4 4 4 2
## [223] 4 2 4 2 4 2 3 4 4 4 2 4 4 2 4 4 4 4 2 4 2 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2
## [260] 4 4 4 4 4 4 1 2 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 2 2 2 2 4 4 4 2 4 2 2 2
## [297] 2 2 2
## Levels: 1 2 3 4

Błąd walidacji krzyżowej jest najniższy przy 4 cechach

randomForest(dane[,-1],as.factor(dane$Kategoria),importance=T)->rf.dane.imp
varImpPlot(rf.dane.imp)

4 najlepsze cechy to V7,V5,V3,V1

randomForest(dane[,c(8,6,4,2)],as.factor(dane$Kategoria))->rf.dane.sel
rf.dane.sel$confusion
##    1  2  3  4 class.error
## 1  7  1 19  3   0.7666667
## 2  1 11  5 49   0.8333333
## 3 10  5 37 14   0.4393939
## 4  1 28 17 91   0.3357664
1-mean(rf.dane.sel$confusion[,5])
## [1] 0.4062099

Używając tej metody osiągamy średnią dokładność równą 42%. Największy błąd jest w klasie 2 i wynosi 82%

##KNN Zacznę od stworzenia danych testowych i treningowych

podział<-createDataPartition(dane[,1],p=0.8,list=F)

dane.train<-(dane[podział,])
dane.test<-(dane[-podział,])

prop.table(table(dane$Kategoria))
## 
##         1         2         3         4 
## 0.1003344 0.2207358 0.2207358 0.4581940
prop.table(table(dane.train$Kategoria))
## 
##         1         2         3         4 
## 0.1083333 0.2125000 0.2208333 0.4583333
prop.table(table(dane.test$Kategoria))
## 
##          1          2          3          4 
## 0.06779661 0.25423729 0.22033898 0.45762712

Teraz zajmę się skalowaniem i standaryzacją

sis_tr<-preProcess(dane.train[,-1])
sis_tr
## Created from 240 samples and 7 variables
## 
## Pre-processing:
##   - centered (7)
##   - ignored (0)
##   - scaled (7)

Żadna cecha nie została pomienięta

dane.train.sc<-predict(sis_tr,dane.train)
apply(dane.train.sc[,-1],2,sd)
## V1 V2 V3 V4 V5 V6 V7 
##  1  1  1  1  1  1  1
summary(dane.train.sc[,-1])
##        V1                 V2                  V3                  V4          
##  Min.   :-1.78861   Min.   :-0.595577   Min.   :-2.014368   Min.   :-2.50869  
##  1st Qu.:-0.77539   1st Qu.:-0.488644   1st Qu.:-0.672388   1st Qu.:-0.53879  
##  Median :-0.09991   Median :-0.344352   Median :-0.001398   Median :-0.02151  
##  Mean   : 0.00000   Mean   : 0.000000   Mean   : 0.000000   Mean   : 0.00000  
##  3rd Qu.: 0.74443   3rd Qu.: 0.009746   3rd Qu.: 0.585718   3rd Qu.: 0.42149  
##  Max.   : 2.85530   Max.   : 7.719756   Max.   : 2.682562   Max.   : 5.09820  
##        V5                V6                 V7         
##  Min.   :-0.3936   Min.   :-5.16093   Min.   :-1.6628  
##  1st Qu.:-0.3936   1st Qu.:-0.57020   1st Qu.:-0.7442  
##  Median :-0.3936   Median : 0.08562   Median :-0.1566  
##  Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.:-0.3936   3rd Qu.: 0.74144   3rd Qu.: 0.9471  
##  Max.   : 7.0154   Max.   : 2.49029   Max.   : 1.8949

Wartość średnia i odchylenie standardowe są równe 0 i 1

set.seed(length(dane.train))
dane.test.sc<-predict(sis_tr,dane.train)
#jako k przjmę pierwiastek z liczby obserwacji
k<-sqrt(length(dane.test.sc$Kategoria))
dane.knn<-knn(dane.train.sc[,-1],dane.test.sc[,-1],dane.train.sc$Kategoria,k)

table(dane.knn,dane.test.sc$Kategoria)
##         
## dane.knn   1   2   3   4
##        1   5   0   3   2
##        2   4   5   1   1
##        3   6   3  20   4
##        4  11  43  29 103
mean(dane.knn==dane.test.sc$Kategoria)
## [1] 0.5541667

Knn nie poradziło sobie najlepiej, średnia dokładność jest róna około 50%.

##Wnioski końcowe (1)Zrobiłem analizę podstawowych statystyk opisowych takich jak wariancja, skośność czy wielomodalność (2)Zanalizowałem macierz korelacji Pearsona, zastosowałem PCA i T-sne do redukcji wymiaru (3)Przeprowadziłęm analizę skupień za pomocą metdy k-średnich, aglomeracyjnej oraz podziałowej (4)Przeprowadziłem kalsyfikację z użyciem lasów losowych (5)Przeprowadziłem klasyfikację za pomocą algorytmu k- najbliższych sąsiadów

Pomimo zastosowania wielu metod statystycznych, analizy danych i klasyfikacji stwierzdam, że wyniki króre otrzymałem pokazują, że modelowanie tych danych może być trudne i mało skuteczne. Żeby osiągnąć zadowalające efekty potrzeba byłoby więcej danych. Trzeba również pamiętać, że używałem danych dotyczących pacjentów z niewydolnością serca, a dane medyczne często nie są możliwe do
jednoznacznej interpretacji.